suppressPackageStartupMessages({
library(DESeq2)
library(dplyr)
library(gplots)
library(ggplot2)
library(here)
library(hyperSpec)
library(limma)
library(readr)
library(tibble)
library(RColorBrewer)
library(VennDiagram)
})
source(here("UPSCb-common/src/R/gopher.R"))
hpal <- colorRampPalette(c("blue","white","red"))(100)
pal <- brewer.pal(8,"Dark2")
suppressPackageStartupMessages(load(here("analysis/batchCorrection/vsd.rda")))
sample.sel <- (vsd$Time %in% paste0("B",3:6) & vsd$Tissue =="SE") |
(vsd$Tissue!="SE" & vsd$Time %in% paste0("B",4:9))
vsd <- vsd[,sample.sel]
vsd$Time <- droplevels(vsd$Time)
# The order is B4-5 ZE; B3 SE; B6 ZE, B4-5 SE, B7-9 ZE, B6 SE
# corresponding to pseudo-time 1-2,3,4,5-6,7-9,10
vsd$PseudoTime <- as.integer(vsd$Time) -1 +
ifelse(vsd$Tissue=="SE",
ifelse(vsd$Time=="B3",3,
ifelse(vsd$Time %in% c("B4","B5"),4,7)),
ifelse(vsd$Time=="B6",1,
ifelse(vsd$Time %in% paste0("B",7:9),3,0)))
goi <- read_tsv(here("doc/MA_list_known_SE_genes.txt"),
col_types=cols(.default=col_character()))
"line_plot" <- function(vsd=vsd,gene_id=gene_id,gene_name=gene_name){
message(paste("Plotting",gene_id))
sel <- grepl(gene_id,rownames(vsd))
stopifnot(sum(sel)==1)
p <- ggplot(bind_cols(as.data.frame(colData(vsd)),
data.frame(value=assay(vsd)[sel,])),
aes(x=PseudoTime,y=value,col=Tissue,group=Tissue)) +
geom_point() + geom_smooth() +
scale_y_continuous(name="VST expression") +
ggtitle(label=paste("Expression for: ",gene_id,"(",gene_name,")"))
suppressMessages(suppressWarnings(plot(p)))
return(NULL)
}
dev.null <- apply(goi,1,function(ro){line_plot(vsd,ro[2],ro[1])})
## Plotting MA_86195g0010
## Plotting MA_98095g0010
## Plotting MA_121578g0010
## Plotting MA_903853g0010
## Plotting MA_52791g0010
## Plotting MA_216909g0010
## Plotting MA_68722g0010
## Plotting MA_147422g0010
## Plotting MA_823242g0010
## Plotting MA_66505g0010
## Plotting MA_66505g0020
## Plotting MA_4929831g0010
## Plotting MA_104203g0010
## Plotting MA_10428962g0010
## Plotting MA_196219g0010
## Plotting MA_8832398g0010
## Plotting MA_42637g0010
## Plotting MA_8990972g0010
## Plotting MA_8227330g0010
## Plotting MA_10429361g0010
## Plotting MA_10434655g0010
## Plotting MA_10432909g0020
## Plotting MA_119006g0010
## Plotting MA_95218g0010
## Plotting MA_10431854g0020
## Plotting MA_255g0010
## Plotting MA_10429533g0010
## Plotting MA_113269g0010
## Plotting MA_166600g0010
## Plotting MA_10163020g0010
## Plotting MA_69685g0020
## Plotting MA_197296g0010
## Plotting MA_58651g0010
## Plotting MA_10435775g0010
## Plotting MA_141759g0010
## Plotting MA_15508g0010
## Plotting MA_9199701g0010
## Plotting MA_39847g0010
## Plotting MA_10428957g0010
## Plotting MA_10436994g0010
## Plotting MA_117784g0010
## Plotting MA_808776g0010
## Plotting MA_8273259g0010
## Plotting MA_10087005g0010
## Plotting MA_10431301g0020
## Plotting MA_176924g0010
## Plotting MA_623831g0010
## Plotting MA_12471g0030
## Plotting MA_102140g0010
## Plotting MA_9725777g0030
## Plotting MA_10429954g0020
## Plotting MA_10437050g0010
## Plotting MA_135287g0010
## Plotting MA_10430507g0030
## Plotting MA_10360330g0010
## Plotting MA_10386237g0010
## Plotting MA_9968141g0010
## Plotting MA_35079g0010
## Plotting MA_1036215g0010
## Plotting MA_10053571g0010
## Plotting MA_1090407g0010
## Plotting MA_17729g0010
## Plotting MA_10434917g0010
## Plotting MA_10437070g0010
## Plotting MA_10435722g0010
## Plotting MA_10435722g0010
## Plotting MA_66505g0010
## Plotting MA_66505g0020
We assume the model to be expression as a function of Stage and Tissue and their interaction i.e. the two variables are not considered independent. However, we will group the Stages, to compare Maturation (SE B3-5 and ZE B3-6) and Desiccation (SE B6 and ZE B7-9).
vsd$Stage<-ifelse(vsd$PseudoTime<7,"Maturation","Dessication")
Stage <- vsd$Stage
Tissue<-vsd$Tissue
design <- model.matrix(~Stage*Tissue)
remove genes with t0o little variance
mat <- assay(vsd)[rowMads(assay(vsd))>0,]
fit <- lmFit(mat, design)
fit <- eBayes(fit)
The possible coefficients
colnames(design)
## [1] "(Intercept)" "StageMaturation"
## [3] "TissueSE" "TissueZE"
## [5] "StageMaturation:TissueSE" "StageMaturation:TissueZE"
fit2 <- contrasts.fit(fit, c(0,0,0,0,-1,1))
fit2 <- eBayes(fit2)
ZSM <- topTable(fit2,p.value=0.01,lfc=0.5,adjust.method="BH",n=nrow(vsd))
sample.sel <- vsd$Stage == "Maturation" & vsd$Tissue != "FMG"
heatmap.2(t(scale(t(mat[rownames(ZSM),sample.sel]))),
distfun=pearson.dist,
hclustfun=function(X){hclust(X,method="ward.D2")},
labRow = NA,trace = "none",
labCol = paste(vsd$Tissue,vsd$Time)[sample.sel],
col=hpal,cexCol=0.8)
fit2 <- contrasts.fit(fit, c(0,0,-1,1,0,0))
fit2 <- eBayes(fit2)
ZSD <- topTable(fit2,p.value=0.01,lfc=0.5,adjust.method="BH",n=nrow(vsd))
sample.sel <- vsd$Stage != "Maturation" & vsd$Tissue != "FMG"
heatmap.2(t(scale(t(mat[rownames(ZSD),sample.sel]))),
distfun=pearson.dist,
hclustfun=function(X){hclust(X,method="ward.D2")},
labRow = NA,trace = "none",
labCol = paste(vsd$Tissue,vsd$Time)[sample.sel],
col=hpal,cexCol=0.8)
grid.newpage()
grid.draw(venn.diagram(list(ZSD=rownames(ZSD),
ZSM=rownames(ZSM)),NULL,fill=pal[1:2]))
GO <- lapply(lapply(list(rownames(ZSD),rownames(ZSM)),
sub,pattern="\\.1$",replacement=""),
gopher,background=sub("\\.1","",rownames(mat)),
task=list("go"),alpha=0.01,
url="pabies")
## Loading required package: tidyr
##
## Attaching package: 'tidyr'
## The following object is masked from 'package:S4Vectors':
##
## expand
## Loading required package: jsonlite
dir.create(here("analysis/limma"),showWarnings=FALSE)
write_delim(ZSD,here("analysis/limma/ZE-vs-SE-desiccation_DE-results.tsv"))
write_delim(ZSM,here("analysis/limma/ZE-vs-SE-maturation_DE-results.tsv"))
write_delim(GO[[1]]$go,here("analysis/limma/ZE-vs-SE-desiccation_DE_GO-results.tsv"))
write_delim(GO[[2]]$go,here("analysis/limma/ZE-vs-SE-maturation_DE_GO-results.tsv"))
write_delim(GO[[1]]$go[,c("id","padj")],here("analysis/limma/ZE-vs-SE-desiccation_DE_GO-for-REVIGO-results.tsv"))
write_delim(GO[[2]]$go[,c("id","padj")],here("analysis/limma/ZE-vs-SE-maturation_DE_GO-for-REVIGO-results.tsv"))
## R version 4.0.3 (2020-10-10)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 18.04.5 LTS
##
## Matrix products: default
## BLAS/LAPACK: /usr/lib/x86_64-linux-gnu/libopenblasp-r0.2.20.so
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] grid parallel stats4 stats graphics grDevices utils
## [8] datasets methods base
##
## other attached packages:
## [1] jsonlite_1.7.1 tidyr_1.1.2
## [3] VennDiagram_1.6.20 futile.logger_1.4.3
## [5] RColorBrewer_1.1-2 tibble_3.0.4
## [7] readr_1.4.0 limma_3.46.0
## [9] hyperSpec_0.99-20201127 xml2_1.3.2
## [11] lattice_0.20-41 here_1.0.0
## [13] ggplot2_3.3.2 gplots_3.1.1
## [15] dplyr_1.0.2 DESeq2_1.30.0
## [17] SummarizedExperiment_1.20.0 Biobase_2.50.0
## [19] MatrixGenerics_1.2.0 matrixStats_0.57.0
## [21] GenomicRanges_1.42.0 GenomeInfoDb_1.26.1
## [23] IRanges_2.24.0 S4Vectors_0.28.0
## [25] BiocGenerics_0.36.0
##
## loaded via a namespace (and not attached):
## [1] nlme_3.1-150 bitops_1.0-6 bit64_4.0.5
## [4] httr_1.4.2 rprojroot_2.0.2 tools_4.0.3
## [7] R6_2.5.0 KernSmooth_2.23-18 DBI_1.1.0
## [10] lazyeval_0.2.2 mgcv_1.8-33 colorspace_2.0-0
## [13] withr_2.3.0 tidyselect_1.1.0 curl_4.3
## [16] bit_4.0.4 compiler_4.0.3 formatR_1.7
## [19] DelayedArray_0.16.0 labeling_0.4.2 caTools_1.18.0
## [22] scales_1.1.1 genefilter_1.72.0 stringr_1.4.0
## [25] digest_0.6.27 rmarkdown_2.5 XVector_0.30.0
## [28] jpeg_0.1-8.1 pkgconfig_2.0.3 htmltools_0.5.0
## [31] highr_0.8 rlang_0.4.9 RSQLite_2.2.1
## [34] farver_2.0.3 generics_0.1.0 BiocParallel_1.24.1
## [37] gtools_3.8.2 RCurl_1.98-1.2 magrittr_2.0.1
## [40] GenomeInfoDbData_1.2.4 Matrix_1.2-18 Rcpp_1.0.5
## [43] munsell_0.5.0 lifecycle_0.2.0 stringi_1.5.3
## [46] yaml_2.2.1 zlibbioc_1.36.0 blob_1.2.1
## [49] crayon_1.3.4 splines_4.0.3 annotate_1.68.0
## [52] hms_0.5.3 locfit_1.5-9.4 knitr_1.30
## [55] pillar_1.4.7 geneplotter_1.68.0 futile.options_1.0.1
## [58] XML_3.99-0.5 glue_1.4.2 evaluate_0.14
## [61] latticeExtra_0.6-29 lambda.r_1.2.4 vctrs_0.3.5
## [64] png_0.1-7 testthat_3.0.0 gtable_0.3.0
## [67] purrr_0.3.4 xfun_0.19 xtable_1.8-4
## [70] survival_3.2-7 AnnotationDbi_1.52.0 memoise_1.1.0
## [73] ellipsis_0.3.1